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Abstract 

We present a robust algorithm for the computation of electromagnetic fields radiated by point sources 
(Hertzian dipoles) in cylindrically stratified media where each layer may exhibit material properties (per¬ 
mittivity, permeability, and conductivity) with uniaxial anisotropy. Analytical expressions are obtained 
based on the spectral representation of the tensor Green’s function based on cylindrical Bessel and Hankel 
eigenfunctions, and extended for layered uniaxial media. Due to the poor scaling of these eigenfunctions for 
extreme arguments and/or orders, direct numerical evaluation of such expressions can produce numerical 
instability, i.e., underflow, overflow, and/or round-off errors under finite precision arithmetic. To circumvent 
these problems, we develop a numerically stable formulation through suitable rescaling of various expres¬ 
sions involved in the computational chain, to yield a robust algorithm for all parameter ranges. Numerical 
results are presented to illustrate the robustness of the formulation including cases of practical interest to 
geophysical exploration. 

Keywords: cylindrically stratified media, anisotropic media. Green’s function, cylindrical coordinates, 
electromagnetic radiation 


1. Introduction 

Analysis of electromagnetic fields in cylindrically stratified media is of great importance in many appli¬ 
cations, such as borehole geophysics mmm- This is a classical problem with separable geometry where the 
components of the tensor Green’s function can be expressed in generic form as 0 Ch. 3],0 



where the integrand factor ^^(p, p') contains various products of cylindrical Bessel and Hankel functions. 
When applicable, such solutions are often preferred to brute-force numerical methods such as finite elements 
and finite difference liiTiisiiiTniiiiiiiiiisiiiiiis] since the former can provide very accurate results with 
computational costs that are orders of magnitude smaller than the latter. This is especially important for 
inverse algorithms relying on repeated forward solutions and which seek to determine sought-after physical 
parameter values (say, layer resistivities) from the knowledge of the field values (measured) at certain 
subterranean locations. 
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However, numerical computations directly based on the canonical expressions of this problem can lead 
to underflow and overflow issues in finite precision arithmetic. This is caused by the poor scaling of cylin¬ 
drical Bessel and Hankel functions for extreme arguments and/or orders, which occur for low frequencies of 
operation and/or extreme values for layer resistivities. In addition, convergence problems in the numerical 
evaluation of the spectral integral on the longitudinal wavenumber kz may occur depending on the sep¬ 
aration distance between the source {p'^(j)'^z') and observation point {p^(j)^z) as well as on the operation 
frequency. To circumvent these problems, a stable formulation based on a suitable analytical conditioning 
of the various factors in the computational chain and a proper choice of deformed integration paths in the 
complex kz plane was recently put forth in [5]. This formulation was shown to be robust to variations on 
physical parameters that span several orders of magnitude. A related formulation to compute static fields 
(electric potentials) due to current electrodes in isotropic layers was described in [16]. 

In this work, we extend the formulation presented in [5] to account for scenarios where the layers compris¬ 
ing the cylindrical stratified media may exhibit anisotropic properties. In borehole geophysics, anisotropy 
is quite common may 

result from geological factors affecting the various Earth layers such as salt water penetrating porous frac¬ 
tured formations and thereby increasing the conductivity in the direction parallel to the fracture and/or 
the presence of clay and sand laminates with directionally dependent resistivities. Here, for generality, we 
assume each layer to be doubly uniaxial, i.e., both the complex permittivity tensor e (which includes the 
conductivity tensor) and the permeability tensor Jl are independently uniaxial, which facilitates the analysis 
of equivalent problems using electromagnetic duality [H Ch. 1]. 


2. Fields in cylindrically-layered uniaxial media 

Most of the basic notation and terminology is adopted from [H Ch. 3]. The section can be regarded as 
a generalization of the formulation presented for isotropic layers in [5| to uniaxial anisotropic layers. 

2.1. General solution in homogeneous, uniaxial media 

Maxwell’s curl equations in uniaxial, homogeneous, and source-free media (with time-harmonic depen¬ 
dence assumed) read as 


V X E = icj/iH, (2) 

V X H = —icjeE, (3) 

where Jl and e are the permeability tensor and complex permittivity tensor, respectively. In the unixial case, 
Jl is written as 


= 


hh 0 
0 Ph 
0 0 


0 

0 

pv 


( 4 ) 


where ph and py are the horizontal and vertical permeabilities, resp. The complex permittivity tensor e 
includes the electric conductivity and it is written as 



0 

o' 



0 

0 

0 


0 

= 

0 


0 

_0 

0 

^v_ 


0 

0 

^p,v T ^ 


where Cp^h and Cp^y are horizontal and vertical permittivities, and ah and ay are horizontal and vertical 
conductivities, resp. In such source-free media, the divergence equations can be written as 


V • (e • E) = 0, 

V • • H) = 0. 
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( 6 ) 

( 7 ) 








Note that in general V • E and V • H in uniaxial and source-free media are nonzero. Indeed, the left hand 
side of § in cylindrical coordinates is written as 




p dp p d(j) dz 
From ^ and (|^, we can obtain 


Similarly, we can obtain 


V-H= 1- 


^v\ 


^h) 

dz J 

^v\ 

dE, 


dz 

Pv^ 

1 dH, 


' dz ■ 




^ • ( 8 ) 


( 9 ) 


To obtain the vector wave equation for E, taking the curl of and using § yields 


V^E- 


1- — 1 V 


.dE, 

dz 


-lUJ 




( 10 ) 


( 11 ) 


where V = Vg + is used and subscript s indicates the transverse components to the ^-component. 
Similarly, the vector wave equation for H can be obtained by taking the curl of ^ and using (10) such that 




l^h 


dz 


= \UJ 


V.s ^ {^h^s + 


( 12 ) 


When the ^-components are extracted from (11) and (12), the equations for the ^-components are written 
as 


- (^1 - ^ + to^fihevE, = 0, (13a) 

^ + Lo^^^ehH, = 0. (13b) 

As usual, Ez and can be solved for using the separation of variables technique. We define the propagation 
constant as /c = oj-s/ph^h with dispersion relation uj‘^ph^h ~ for the longitudinal (or vertical) kz and 

transverse (or radial) kp wavenumbers. Two different anisotropic ratios can be defined in such media: the 

anisotropy ratio for the complex permittivity as and the anisotropy ratio for permeability as 

= \ —• It is also convenient to define two scaled radial wavenumbers as = — and ko = —. With 
the above definitions, the general solution to the vector wave equation in such media becomes 


Ez 

Hz 


AnJn (fcpp) + (fcp/o) 

CnJn (kpp) + (kpp) 


^inp^ikzZ ^ 


(14a) 

(14b) 


with An, Cn, and determined by boundary conditions. The dispersion relations for kp and kp are 

— {io^Ph<^h - kl) = kp, (15a) 

— {J^ph^h - kl) = Pp. (15b) 
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Figure 1; Two different cases of two uniaxial cylindrical layers with relevant reflection and transmission 
coefficients in the /9^-plane; (a) Outgoing-wave case and (b) Standing-wave case. 


The transverse {p and <p) field components can be expressed in terms of the above longitudinal components 
directly by using Maxwell’s equations |1], to yield 


E. = 

H, 


1 


[ dz 

1 


^ dE, . . ^ „ 

Vs— - lojphZ X VsHz 


or, in a convenient matrix form. 


OH, . ^ ^ 

Vs—-h lujehz X VgEz 

oz 


= p [\kzVsEz - iujphZ X VsHz], 

= p [ik^VsH^ -F iuje^z x VgE^]. 


(16a) 

(16b) 
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(17a) 

(17b) 


It should be noted that and depend on horizontal medium properties and fih- 


2.2. Local reflection and transmission coefficients 

When a number of cylindrical layers with different properties are present, the appropriate boundary 
conditions need to be enforced into the solutions. This is typically done via reflections and transmissions 
coefficients. In this section, local reflection and transmission coefficients are first derived for the two-layer 
case. Later, this will be extended to the case with arbitrary number of layers. 

The local coefficients can be classified into two types: out going-wave type and standing-wave type 
depending on the relative location of the source versus the observation point, as illustrated in Fig. As the 
general solution of Ez and Hz for uniaxial media are slightly different from those for isotropic media, so are 
the local reflection and transmission coefficients. Nevertheless, the expressions for the generalized reflection 
and transmission coefficients (to account for more than two layers) in terms of local ones remain the same 
as those in [5]. 
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2.2.1. Outgoing-wave case 


Based on ( |14a[ ) and ( |14b[ ), outgoing waves in a uniaxial medium can be expressed as 

= H W(fcpp) ■ a, 


'eJ 


'H^n\kpp) 0 

ez 



0 H!^l\kpp)_ 

hz_ 


( 18 ) 


where the column vector a includes dependence. Since the source is embedded in layer 1 for the 

outgoing-wave case depicted in Fig. la, the ^-components of the total fields in layer 1 and layer 2 is expressed 
as 


E. 


' = H W (fcipp) ■ ai + J (fcip/9) ■ R12 ■ ai, 

= ■ Ti 2 ■ ai. 


■z2 


(19a) 

(19b) 


Likewise, using (17b), we have 


E^i 

H^2 

-£’<^2 


“ '^1 ^ 4>n{klpP) •R-12 '^1, 


= H2(£2p/9)-Ti2-ai. 


(20a) 

(20b) 


From ( |19a[ ) through ( |20b[ ), two types of matrices only depending on p are defined as 


^ zniJ^ipP) — 


B0n(^ipP) — j^2 


1 


0 

0 ^niJ^ipP) 

iujehikippB'^{kipp) -nk^Bnikipp) 
nk^B^ii^kipp^ ’\uj PhjkippB^ikipP^ 


(21a) 

(21b) 


where B^ is either Btn"^ or kip = uj‘^phEhi — and e^i and phi a-re the horizontal complex permittivity 
and permeability in layer respectively. Applying the pertinent boundary conditions, viz., continuity of 2 (- 
and ^-components at p = ai, to di^-d^ yields 

H[b(fcipai) + J2„(fcipai) ■ Ri2 • ai = H[j,^(fc2pai) • T 12 ■ ai. 


H^^i(fcipai) + J0n(fcipai) • R 12 I • ai = H^^^(fc2pai) • T 12 • ai. 
For notational simplicity, a shorthand notation is defined such that 

^ani^ip^j) — ^aij‘ 


(1)/ 


(22a) 
(22b) 

(23) 

In the right hand side of ( [^ , the first, second, and third subscripts indicate the relevant components of the 
fields, radial wavenumbers, and radial distances, respectively. For notational simplicity, the Hankel function 
superscript and subscript (kind and modal number) are suppressed in the following. Consequently, (22a) 
and (22b) are re-expressed as 


[^zii E Jzii • R12] = iiz2i • T 12, 
[H^ii H- • R12] = H021 • T12. 
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From (24a) and (24b), we obtain 



_ _^ _ 

-1 

r 

Ri2 = 

J;2ll — ^z21 ' H021 • J011 


H.21-H^21 

HI 

to 

II 

^z21 — Jzll • J^ll • H021 

-1 

^zll — Jzll 


• — tLzii 


(25a) 

(25b) 


2.2.2. Standing-wave case 

For the standing-wave case depicted in Fig. we now have 


and, using (17b), 


Ezi 

Hzi 

Ez2 

Hz2 


E(f)i 

H^2 

Ecf)2 


J zn {kipp) • T 21 • a2, 

H 2^1 {k2pp) • R 21 • a 2 + J zn{k2pp) ■ a2, 


J 0n(^lpP) ■ T 21 ' a2. 


H^^2(^2p/3) • R 21 • a 2 + J 4 m{k 2 pp) ■ a 2 . 


Applying the boundary conditions at p = ai to yields 


HW(fc2pai)-R 21 + J zn{k2pai) 

( 1 ) 

(pn 


J zn{kipai) • T21 • a2 — 

J pnikipai) • T 21 • a2 = H ^}^{k2pai) • R 21 + J 0n(^2p<^i) 


• a2, 

• a2, 


(26a) 

(26b) 

(27a) 

(27b) 

(28a) 

(28b) 


which can be re-expressed as 


Jzll • T 21 — [H;221 ' 1^21 + J;22i] 5 

J^ll • T 21 = [H021 -1^21+ J(/)2l] • 

Therefore, R 21 and T 21 are written as 


R 21 — 
T 21 = 


H,21 

Jzll - 


■ • J(/)ll ’ H 


021 


tiz21 • H021 • J011 


-1 


-1 


' J^ll ’ J021 


fz21 


»z21 


■ tiz21 • H021 • J021 


(29a) 

(29b) 


(30a) 

(30b) 


Furthermore, the local reflection and transmission coefficients (25a), (25b), (30a), and (30b) can be succinctly 
rewritten as 


to 

II 

01 

^ J- 

^z21 ’ H021 • — iizll 

_ _^ 

_^ - 


R 21 = • 

J;2ll • J^ll • J021 — Jz21 

5 

ti 2=Dy ■ 

J^ll ' J^ll ' H^ll 

HI 

to 

II 

01 

^ J- 

J;221 ~ H^21 • H021 • J021 


(31a) 

(31b) 

(31c) 

(31d) 


6 


















































where 


(32a) 

(32b) 


II 

IQ 

— H;s21 • H021 • 


- _ _ _ ^ _ - 

Db = 

H;221 — ^z\\ ' • H021 


2.3. Spectral representation of the Greenes function 

In this section, we derive convenient analytical expressions for the Green’s function expressed in cylindri¬ 
cal coordinates that will facilitate further analysis in such uniaxial media. Specifically, we obtain expressions 
for the ^-components of the fields produced by a point source (arbitrarily-oriented Hertzian electric dipole). 
In this case, the source writes as 


J(r) = Ila 5 {y — r'), 

where II is the dipole moment. By taking the curl of Faraday’s law, we obtain 

V X V X E = iuV X pH. 


With the presence of the source, (34) becomes 


V^E- ( 1 - —) = -iwV X ^H +V 

' eh J az \ lioeu 


(33) 

(34) 

(35) 


where the divergence of E is deduced from ^ and the continuity equation V • J — iuupv = 0 is applied. 
Extracting the 2 )-component of the above, we can easily show that 


'V Ez uj phCyEz — ( 1 - - 


ChJ dz‘^ 


d^Ez . . , 5 /V-J 

= -lUOPhZ • J + — T 


Using (33) and e^^^^ dependence shown in (14a), we obtain 


V^E, + k‘^E, = - — 


d 


oz 


dz \ iuJCh 


(5(r - r'), 


where 


= J^Hhev + ( 1 - -]kl = — {uj'^liheh - kl) +kl = kl + kl, 


7 2 2 
k =UJ PhCh- 

Therefore, Ez is obtained via the scalar Green’s function. 


E. = 


ill 

U)€h 


d 


k^iz-a') + —V 


where 


g{r - r') = 


pi/c|r-r'| 


^(r - r'), 


(36) 

(37) 

(38) 

(39) 

(40) 

(41) 


47r|r — r'l ’ 

The derivation for the 2 )-component of the magnetic field can be done quite similarly by taking the curl 


of Ampere’s law. The equivalent to (36) becomes 


V^HzEuo^p^ehHz-il-^ 






lih J 


-5 • V X J. 


(42) 
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Since Hz also has dependence as shown in (14b), (42) reduces to 


where 


+ Ph^ = -Ilz • V X a'5{v - r'), 


P = + [l - — ] P^ = — {Pnheh - Pz) +Pz=Pp + Pz- 


l^h 


l^h 


Now, Hz is obtained via 


where 


Hz = -Ilz ■ V' X a'g{r - r'), 

gi/c|r-r'| 


g{r - r') = 


47r|r — r'l 


(43) 

(44) 

(45) 

(46) 


The 2 )-coniponents of electromagnetic fields expressed in (40) and (45) can be expanded as the linear 
combination of all spectral components. In other words, using the spectral representation of the scalar 
Green’s function. 


gi/c|r-r'| ((f)-(f)') poo 

|r-r'| ^ ^ 2 

' ' n= —oo 


/ oo 

-OO 


(47) 


Ez and Hz in homogeneous, uniaxial media can be written as 


Ez 

Hz 


HI 


IlTLOeh 


OO ^OO 

n— ^ J —(X) 


Jn{kpP<)Hfa\kpP^) 


0 Jn{kpp<)H^^\kpp>) 


where h' is an operator acting on the primed variables on the left and written as 


2 


izP + £.V)-a' 

iuj€h^'' z X 'V' 


(48) 


(49) 


In order to account for multilayers, let us consider cylindrically stratified media with the source in layer 
j. Then, the reflection terms from the boundaries at p = aj and p = Uj-i are added to (48) such that 


E: 


^3 


Hz 


where 


HI 


lTTU)€hj 


^ POO , _ 

^ {Jzjp, ■ nzjp, + Uzjp ■ a,„(pO + J,,-, ■ b,n(pO} ■ 5', 

n— ^ J —(X) 


^zjpc ' ^zjp> — 


^zjp — 


^zjp — 




Jnikjpppufi ikjpp>) 0 

0 Jn{kjpP^)Hn {kjpPy^ 

Hff\kjpp) 0 

0 Hff\kjpp) 

Jn{kjpp) 0 

0 Jni'kjpP') 

{zk] - ikzV') ■ a' 
iuiehja' • z X 'V' 


(50) 

(51a) 

(51b) 

(51c) 

(51d) 
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Recall that notations of (51a), (51b), and (51c) are based on (21a) and (23). In (51d), kj = and 

jj^hj and Chj represent horizontal permeability and complex permittivity in layer j, respectively. When a' is 
represented in cylindrical coordinates as a' = p'apf + + z'az ', 






1 

1 

+ 

’ ' 

d 

T —— 

ikzOLp' 

1 

0 


nujehj 

P' ^ . 

dp' 

\UJ€.hj OL(f)' 


(52) 


Using the two constraint conditions at p = aj and p = aj_i detailed in [4], two unknowns and 

hjn{p') appeared in (50) can be determined such that 


II 

^ Ri,i-i ■ Ri,i+i 

’ R^j—1 

^zjp' + R jy+i • ^zjp' 

bjn — 

^ ~ R j,i+i ' Ri,i-i 

' Rj^'+i ' 

^ zjp' + R • ^zjp' 


(53a) 

(53b) 


Note that the second brackets in the right hand sides of (53a) and (53b) are slightly different from those for 
isotropic media. When p > p' ^ p^ = p' and p> = p. The curly bracket in (50) is expressed as 




H 


ZJP 


'^zjp T J 


^jn{p ) T ^zjp ’ t)jn(p ) 


zjp ' 


•M,+ . 


^zjp' + • H 


zjp' 


(54) 


On the other hand, when p < p' ^ p^ = p and p^ = pk The curly bracket (50) is now expressed as 


^zjp< ’ ^zjp ' ^jn{P ) T ^ zjp ’ t>jfn(P ) 


^zjp + ^zjp ' 


M,_ 


+ R 




(55) 


Again, (54) and (55) are slightly different from those for isotropic media. When the field layer is not the 
same as the source layer, the approach is the same as that in [4]. In summary. 


E, 

Hz 


ill 


Airuehj ^ 

77 ,= — 


^ pOO , 

i J-oo 


(56) 


where: 

for Case 1: p and p' are in the same region and p > p'. 


EniP^ P ) 


'^zjp T J zjp ’ 


zjp ' ^j,j-\-l 


■Mj+- 


Jzjp' + ' H 


^zjp' 


(57a) 


for Case 2: p and p' are in the same region and p < p'. 


Efiip^ p) — 


•zjp 


^zjp ' R j,j-l 




H 


ZJP' 


R 


j,j + l • ^zjp' 


(57b) 


for Case 3: p and p' are in different regions and p > p'. 


Efiip^ P ) — 


^^zip T ^z 


zip ’ R i,i+l 


•N 




M 


i+ 


^zjp' + R ' ^zjp' 


(57c) 


for Case 4: p and p' are in different regions and p < p'. 


Enip^ P ) — 


■ ^^zip ' R i,i — l 




^zjp' + R j,j+l ■ ^zjp' 


(57d) 
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Table 1: Definition of range-conditioned cylindrical functions for uniaxial complex permittivity media. 


Small arguments 


Moderate arguments 


Large arguments 


Jn^kipCLj^ — Gi(ljJji{kip(lj^ 

Hk^\kipaj) = G-^apH^^\kipaj) 
Hn^^\kipaj) = Gr^aJ^Hn^^\kipaj) 


— Pij JniJ^ipCl'j^ 

Jn{kipaj) = PijJ^{kipaj) 


H^^\kipaj) = Pr^M^\kipaj) 

Hn^^\kipaj) = 


Jnikip^j) — ^ ^ Jni^kipdj^ 

Jnikipaj) = e'^^”Mu'„{kipaP 
(kipaj) = H^rl^ {hpaj ) 

Hn^^\kipaj) = e-'^>^H'Ji^^(kipaP 


3. Range-conditioned formulation 

As noted before, the poor scaling behavior of Bessel and Hankel functions for extreme arguments and/or 
orders causes instabilities in the numerical computation of the field expressions above under some param¬ 
eter ranges. This section discusses how to stabilize the computation and provides relevant mathematical 
derivations. 


3.1. Range-conditioned cylindrical functions and matrices 

Range-conditioned cylindrical functions derived in [5] are modified for uniaxial media because two dif¬ 
ferent scaled radial wavenumbers kip and kip appear as part of the function arguments. Here, the subscript 
i is the layer index, and, to recall, Ki^ = yjcfffif^i and Kip = a/ Rm/R vi are the anisotropy ratios of com¬ 
plex permittivity and permeability in layer i, respectively. Here, we will only focus on those aspects that 
differ from [5] . The reader is referred to [5] for the fundamentals of this stabilization approach. Table 
shows the definitions of the range-conditioned cylindrical functions, indicated by a hat, for uniaxial complex 
permittivity media, where: 





(58) 


P 

-L tn 


1, ^ \i\Jn{kipaj)\ ^ <Tr 

\Jn{kip0^p\i if \ Jn{kipO'p\ ^ -^r 


(59) 


kip = '^rn 


ki 


ip 


k'ip + i^"p 


{k'ip + ik'lf - i/c" ) 

+ (<e) 


!<.' h" _ k" y 

^ie^ip ^ie^ip 


(60) 


where is the magnitude threshold for moderate arguments [5]. Note that subscripts i and j are arbitrary. 
Range-conditioned functions for uniaxial permeability media can be similarly constructed. The multiplica¬ 
tive factors associated with the new functions shown in Table can be classified into two types: a-type and 
/3-type, whereby the relationship between the original cylindrical functions and the range-conditioned ones 
can be succinctly expressed as 


Jn{kip^j^ — ^ij Jnikip^j^ •) 

(61a) 


(61b) 

(kipaj) = aij (kipaj) , 

(61c) 

H'Ji^\kipaP = aijH'Ji^\kipap. 

(61d) 
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Table 2: Definition of aij and I3ij. 


Argument type 

aij 

flij 

Small 

GPaA 


Moderate 


Pij 

Large 




The definitions of otij and ^ij are provided in Table Similarly to the isotropic case, a^j and ^ij exhibit 
two important properties to ensure a stable computation [5]: 

1. Reciprocity. 

au = (62) 


2. Boundness. 




, eXj, 


< 1, for Om < dn 


(63) 


For the anisotropic case, it is convenient to also derive range-conditioned cylindrical ‘matrices’ because 
2x2 matrices rather than scalar factors appear in the computation of the reflection and transmission coef¬ 
ficients in layered media. We use hats to denote those matrices as well. From (21a) and (21b), we obtain 




Jnifip^j ) 

1 

o 


flij 

1 

o 

1 

O 

Jniflip^j ) 


0 

flij 

H^n\hpaj) 

0 



1 

1 

o 


,aj) 


0 


— Pij ‘ ^ zn — ^ zn ’ f^ij •> 


0 

O^ij 


— a - • • H — H • a - • 


J (pnikipaj) — ^2 


K'nikipa,) = 


ip^j 


kfpaj 


^^pUjjU nyr^ip'^j ) 

^zJn(J^ip(^j ) 




nkzJfiiflip^j ) 



0 

iuj llhi’kip^^j ) 


0 

flij 

nk^^Hfi {kipUj 

) 


-1 

p^ 


= Jc 


-nkzH^^\kipaj) -iLOiJ,hi'kipajH'J:^\kipaj) 


0 a 


V 


' ^ii ’ 


(64a) 

(64b) 

(64c) 


= Hg-a,,. (64d) 


Note that two matrices in (64a) and (|64b|) are diagonal, so they commute. 


3.2. Range-conditioned reflection and transmission coefficients 

Using the redefined matrices in (64a) - (64d), local reflection and transmission coefficients for the two- 
layer medium depicted in Fig. [^can be also redefined. First of all, two intermediate matrices. Da and D^, 
are redefined as 



^ —1 ^ 


A ^ ^-1 ^ 


II 

IQ 

Jzll • — ^z21 • 0^21 • 0^21 ' ^021 ' • f^ll 

= 

Jzll — ^z21 • H021 • 

’ ^11 — • fJll^ 





(65a) 

Db = 

D.z21 • 0:21 — Jzll • • H021 • <^21 

= 

^z21 — Jzll • J(f)ll • H021 

• ^21 = Db • ^21- 


(65b) 
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Figure 2: Reflection and transmission coefficients for two cylindrical layers. 


Therefore, the local reflection and transmission coefficient are redefined as 


R12 — D /I • 



_ 1 ^-1 

^ —1 ^ 

H;221 • H(/)21 ’ Hc/)11 — H;2ii 

= Pll ■ ■ 

^z2i • <^21 • OL 21 • H021 • H^ii • an — • an 


^-1 

ail • Da 


^-1 ^ ^ 
Hz21 • H021 • H^ii — Hzll 


■ OLii — ail • Ri2 • 0^11, 


(66a) 


R 21 — • 


Jzll • J^ll • J021 — Jz21 


-1 


= ^21 • • 


^zll ' ' J^ll ' J(/)21 • [^21 ^z21 ' ^21 


— P 21 ’ 


Jzll • J(/,ll • J(/)21 — Jz21 


' P 2 I — P 2 I ' R 2 I • ^21^ 


(66b) 


ti 2=Dy 




^ ^ —1 ^ 

— J;^!! • J^n • H^n 

= OL21 ■ T>b ■ 

^zii • ^11 — Jzii • ^11 • • J^n • H^n • an 


— ^21 ’ 


H;2ll “ • J(/)ll • H^II 


• <^11 — ^21 ’ T 12 • 0 : 11 , 


(66c) 



_ 1 ^-1 

^ ^ ^-1 ^ _ 

Jz21 — H;221 • H(/)21 ' J021 

= /3ll ■ ■ 

J;221 • ^21 ~ ^z21 • <^21 * <^21 ’ ^^21 ' J(/)21 ' ^21 


^-1 

= ail • 


J .21 - H 


z21 • H021 ' 




• f^21 — <^11 • T 21 • /321. 


(66d) 


We can proceed to redefine generalized reflection coefficients for multilayers, which are functions of local 
reflection and transmission coefficients. The generalized reflection coefficient for the out going-wave case for 
three cylindrical layers depicted in Fig. [^is modified to 


Ri 2 = Ri 2 + T 21 • R 23 • — R 2 I • R 23 ] • T 12 


— ail • |R 12 + T 21 • ^21 ’ ^22 • R 23 • ^21 ’ <^22 • I — R 21 • ^21 ’ <^22 • R 23 • ^21 ' ^22 

= ail • R 12 • OLii. 


-tiaj 


an 

(67) 


Note that the magnitude of the multiplicative factor /32i 


a 22 in (67) is never greater than one due to 


the boundness property, which guarantees moderate magnitude for R 12 in any case. Since the associated 
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Figure 3: Generalized reflection coefficients for three cylindrical layers: (a) R 12 for the outgoing-wave case 
and (b) R 32 for the standing-wave case. 


multiplicative factors an shown in (67) are the same as those for R 12 (see ( | 66 a| )), they do not change when 
more than three layers are present. Therefore, the redefined generalized reflection coefficient between two 
arbitrarily-indexed adjacent layers for the out going-wave case can be expressed in general as 




( 68 ) 


where 




-1 


2 , 2+1 • 


(69) 


The generalized reflection coefficient for the standing-wave case for three cylindrical layers depicted in Fig. 
I3blis modified to 


R 32 = R 32 

= /^32 


■{ 


• T 23 • R 21 • [l — R 23 • R 21 ] ' T 32 

R 32 + T 23 • ^21 ' ^22 * R 21 * ^21 ' ^22 * I “ R 23 ' ^21 ' ^22 * ^ 21 ' ^21 ' ^22 


— f^32 '1^32 • 


32- 


• T 32 I • [3^2 

(70) 


Again, the multiplicative factor / 32 i •a 22 in (70) is never greater than one in magnitude due to the boundness 
property, which also guarantees moderate magnitude for R 32 in any case. Again, when more than three 


layers are present, the multiplicative factors in (70) are not affected. Therefore, the redefined generalized 


reflection coefficient between two arbitrarily-indexed adjacent layers for the standing-wave case is expressed 


as 


where 


1 ^ 2 + 1 ,2 /^ 2 + l ,2 ’ 1^2 + 1,2 ' /^ 2 + 1,25 


1 ^ 2 + 1,2 1 ^ 2 + 1,2 ~b T 2 , 2+1 ’ /^ 2 , 2—1 ’ ^22 ' ^^ 2 , 2—1 ’ /^ 2,2 —1 ’ ^22 

I I^2,2 + l ’ /^2,2—1 ’ " ^2,2 — 1 ' /^2,2—1 ’ 


-1 


i + l ,+ 


(71) 


(72) 
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Figure 4: 5'-coefficients for three cylindrical layers: (a) S 12 for the outgoing-wave case and (b) S 32 for the 
standing-wave case. 

Before we proceed to obtain generalized transmission coefficients, the redefinition of the so-called S- 
coefhcients [H Ch. 3] is necessary, as they represent local transmission factors in the presence of multilayers. 
The ^'-coefficient for the out going-wave case for three cylindrical layers depicted in Fig. I^is modified to 


S12 — [l — R 21 • R 23] • T12 


= /3: 


21 


1 -1 


T 12 • ail 


I — R 21 • P 21 ' ^22 • R 23 • P 21 ' ^22 
= ^21 ' S 12 • ail. 

Therefore, the redefined arbitrarily-indexed ^'-coefficient for the outgoing-wave case is written as 

S 2 , 2+1 ~ /^ 2 + l ,2 ’ ^ 2 , 2+1 ' ^225 


where 


i,i+l — 


I — R 2+1,2 • /3i + l,i • <^2 + 1,2+1 • R 2+1,2 + 2 • /^i + l,^ • <3:^ + 1,i + 1 


-1 


2 , 2+1 • 


(73) 

(74) 

(75) 


The ^'-coefficient for the standing-wave case for three cylindrical layers depicted in Fig. 4b is modified 


to 


S 32 = [l — R 23 • R 21 ] • T 


32 


— a22 


1 -1 


' T 32 • /3; 


32 


I — R 23 • ^21 ' ^22 • R 2 I • ^21 ' ^22 

= a22 • S 32 • 

As a result, the redefined arbitrarily-indexed ^'-coefficient for the standing-wave case is written as 

^ 2+1,2 ^22 ' ^ 2 + 1,2 ’ /^2+1,25 

where 


S 2+1,2 — 


I r^2,2 + l ' /^2,2— 1 ’ ' ^2,2 — 1 ' /^2,2—1 ’ 


■ 2 + 1,2- 


(76) 

(77) 

(78) 


Let us now consider the generalized transmission coefficient for the outgoing-wave case {i > j) in cylindrically 
stratified media, which is expressed as 


T — T • S i-2,2-1 




(79) 
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(79) can be modified in a way that 


ji — T i — \^i • S 2—2,i—1 • • • S 

~ ^ i-l,i ’ I ’ ^/c+l,/c+l ’ ^ k,k-\-l | ’ OLj 


\k=j 


~ f^iA-1 ji ’ ^jj' 


(80) 

(81) 


The magnitude of the multiplicative factors • cxk-\-i^k-\-i in (80) is never greater than one, which 


stabilizes the computation of T ji. The product in ( |8Q| ) is the product of a number of 2x2 matrices, so the 
order of the product should be specified. The 2x2 matrix for k = j and 2x2 matrix for k = i — 2 should 
be placed in the rightmost and leftmost in the matrix product, respectively. Furthermore, when i = j + 1, 
the matrix product reduces to an identity matrix. It should be also noted that the associated multiplicative 


factors shown in (81) are the generalized version of those shown in (66c). 


Next, the generalized transmission coefficient for the standing-wave case {i < j) in cylindrically stratified 
media is expressed as 


Tii = T 


2 + 1,2 


■S 


2 + 2 , 2+1 


''' ^id-i* 


(82) 


Similarly, (82) is modified to 


= T, 


2 + 1,2 ' S ^ + 2,2 + 1 • • • S jj-l 

^ / i-1 _ ^ ^ 

= OLii • T ’ I f^k,k — l ’ ^kk ’ ^ k-\-l,k 

\/c=2+l y 


1 


— ' T jf2 • f3j j_-^. 


(83) 

(84) 


Again, the magnitudes of the multiplicative factors • ^kk in (83) are never greater than one. For the 

matrix product in (83), the 2x2 matrix for /c = i + 1 and 2x2 matrix for k = j — 1 should be placed in 
the leftmost and rightmost, which is opposite to the outgoing-wave case. Furthermore, when j = i -h 1, the 
matrix product reduces to an identity matrix. The associated multiplicative factors shown in (84) are the 
generalized version of those shown in (66d). 


Several auxiliary coefficients appeared in (57a) - (57d) should be redefined properly as well. For the first 


integrand type, shown in (57a), M is redefined as 


— 

r_ ^ ^ 1 

-1 

_ A _ A 

M,+ = 

I ~ j,i+i 

= 

^ ~ ’ ^+i-l ’ ^jj-l ■ ^33 ■ ^ij + l ’ ^33 




n -1 


^ ^i,i-i ’ ^i4i-i»j] ' ^id-i ' ' ^jj ' ^i,i+i ' ^i,[i-i,j] ' ^jj 






(85) 


where the radial distance corresponding to subscript [j — 1, j] is = caj-i + (1 — c)aj, 0 < c < 1. Two 

extreme choices of = aj-i and = aj) can be used for notational convenience but 

these are not useful in the redefinition of the integrand, as clarified below in Section [3^ 
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For the second integrand type, shown in (57b), M is redefined as 


— 

r_ ^ ^ 1 

-1 

A _ A _ ' 

II 

1 

I “ j,i+i ■ 

= 

^ ~ ^33 ■ + l ’ ^33 ’ 


-i -1 


= OL 






-1 






( 86 ) 


Again, the two extreme cases of are undesired for the proper redefinition of the integrand as shown 

in Section [3^ 


For the third integrand type, shown in (57c), N is redefined as 



r_ _ ^ 1 

-1 

^ _ A 

Ni+ = 

I 1 ■ 

= 

I Pi,i — 1 ’ — 1 ' Pip—1 ’ ^ii ’ ’ ^ii 


n -1 




1 -1 


I 1 ' Pi,i — 1 ' ^ii ' ' f^i^i — 1 ’ ^ii 


■ OLi^i-i 


■ •Ni+ 


(87) 


Finally, for the fourth integrand type, shown in (57d), N is redefined as 

— n -1 


N,_ = 


= CXi 


I OLii • ' ^ii ’ f^i^i — 1 ’ —1 ’ f^i^i — 1 

-1 




I 1 

I ‘ 1^2,2 —1 ’ 1 ’ ^ 

= aii •Ni_ ^ 

3.3. Range-conditioned integrand 

For Case 1 in (57a), there are four arguments of interest: kjpOj-i^ kjpp'^ kjpp, and kjpOj. For convenience, 
we let aj-i = ai, p' = a 2 , p = as, and Oj = a 4 so that ai < a 2 < as < a 4 . The integrand is redefined as 


( 88 ) 


Fn (/>,/>') = 


^zjp + ^zjp ' R j,j + l 


M 


3 + 


^zjp' + • ^zjp' 

• ^j2 • M • a^-2 


^zjp ' ^j3 + ^zjp ' f^j3 ’ ^i4 • R j,j+l • OLj4 

Jzjp' • ^j2 + ^jl ' ' ^jl ’ ^zjp' • OLj2 

j2 ’ ^j3 ' ^zjp T f^j3 ' ^i4 ' ^zjp ' ' Pj2 ’ ^i4 

(3jl • a.j2 ’ Hjj-I • Hzjp' • Pjl • ^j2 


(89) 


M 


3 + 


^zjp' 


(90) 


Note that, as (89) shows, the corresponding radial distance to the subscript [j — l,jf] in (85) is chosen to 
be a 2 , neither ai nor a 4 . To be more specific, the radial distance of the source p' is selected. The choice 
enables the left and right squared bracket factors in (90) to be balanced and yields a stable computation. 
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For Case 2 in (57b), four arguments are of interest: kjpaj-i^ kjpp, kjpp'^ and kjpaj. Similarly, we let 
aj-i = ai, p = a 2 , p' = < 23 , and aj = so that ai < a 2 < as < a^. The integrand is redefined as 




^zjp + ^zjp • 






' ^zjp' 


Jzjp • ^j2 T ^Zjp ' OLj2 • /3ji • • /3 




• a^-3 • M^_ -/^ 


H 




Cljfs + ^zjp' ’ f^jS 


^j2 ' • Jzjp + ^jl • <^j2 • H2;jp • Hjj-i • • CXjS 

^zjp' + f^js ' ^j4 • • f^js ' ^j4 


i3 




(91) 


(92) 


It should be noted that, as (91) shows, the corresponding radial distance to the subscript [j — l,jf] in ( 86 ) 
is chosen to be a 3 , the radial distance of the source (neither ai nor a^). Again, this choice enables the left 
and right squared bracket factors in (92) to be balanced and yields a stable computation. 

For Case 3 in (57c), there are 6 arguments of interest: kipai-i^ kipp, kipai^ kjpaj-i^ kjpp', and kjpaj. 
We let a^_i = ai, p = a 2 , = a 3 , aj-i = 5i, p' = 62 , and aj = 63 so that ai < a 2 < a 3 and 61 < 62 < 63 . 

The integrand is redefined as 


Fn(P: P ) 


• N • T • M • 


^zjp' T R ?,?—1 ■ H, 


^zip T ^zip ’ Rzji+l 

^zip ’ ^i2 T ^zip (^i2 ' ’ Ri,2+1 ’ 

■ (/3ii ■ N i+ ■ an 

^zjp' ' f^j2 T f^jl ’ Ri,i-1 ' f^jl ’ '^zjp' ’ OLj2 


ij'-i ■ ^zjp' 


) ' il^il ji ’ ’ (^j2 ' ^ j+ ' 


(93) 


Pil ' ^i2 ' ^zip T Pi2 ' ^iS ' ^zip ' Ri,i+1 ' Pn ' ^i3 


• • T • f3j2 ' ^jS ' M j + 


^zjp' T f^jl ’ ^i2 ’ R j,j —1 ’ ^zjp' ' f^jl ’ ^i2 


(94) 


In should be stressed that the corresponding radial distance for M in (93) is now 62 , which is the radial 
distance of the source. 

For Case 4 in (57d), the arguments of interest are the same as those for Case 3. The integrand is redefined 
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as 


Fn(P5 p) — 


^zip “1“ '^zip ' 1 ■ N 2— • T jf2 ’ ' '^zjp' H“ R'jfjjf+l ’ ^zjp' 

^zip ' Pi2 '^zip ’ ^i2 ' f^il ’ 1 ' f^n 


(ai3 • Ni_ • • [ ai3 • T ji • • M 


^^zjp' ' ^j2 “}“ ' ^i3 ■ ^zjp' ' j2 


( 95 ) 


(^i2 ' ^*3 ■ ^zip “1“ Pil ' ^i2 ' ^^zip ' 1 ' Pn ' ^i3 


■ N • T • Pji ' ^j2 ' j — 


'^zjp' + /3j 2 ' ^i3 ■ ■ ^zjp' ■ /3j 2 ' ^i3 


(96) 


Again, the radial distance of the source &2 is chosen for the corresponding radial distance for Min (95) 


3.4-. Azimuth modal summation 

The spectral representations of electromagnetic fields involve an infinite series as the azimuthal summa¬ 
tion. The three components of the electromagnetic fields are expressed as [S] 


Ez 

H, 

Ep 

Hp 

E^, 

Hj, 


ill 


4:TTU€m J- 


:/ 

1 J-C 




OO 

^ ein(^-^')F„(p,p')-6' 


ill 


/ 


Al. pikz(z-z') }_ 


OO 

^ ■ L„(p) ■ M„ ■ R„(pO ■ 6; 


ill 


■iTTUChj J-oc 


/ OO -j ^ , 

dk,pkz(z-z') ■ L„(p) ■ M„ ■ R„(p') ■ 5 ' 

-OO 


(97a 


(97b) 


(97c) 


To expedite the computation, it is possible to fold the series above by exploiting symmetries of the cylindrical 
eigenfunctions so that only zero and positive orders remain. The expressions of folded sums are quite similar 
to those for isotropic media in [5], so final results are not provided here. In addition, the spectral integrals 
(97a), (97b), and (97c) typically cannot be computed along the real axis in a robust fashion. These integrals 


should instead be numerically evaluated along suitably deformed integration paths in the complex kz plane. 
To this end, either the so-called Sommerfeld integration path (SIP) and the deformed SIP (DSIP) can be 
used with the optimal choice among the two, depending on the longitudinal distance between source and 
observation points [5]. Moreover, it should be noted that the concept of the direct field subtraction when 
source and observation points are in the same layer, discussed in [5], still applies to uniaxial media. As 
analytical field expressions in such media are available only in certain cases (see Appendix A), direct field 
terms for isotropic media are used in the computation. 


4. Numerical validation results 


This section provides some validation results for the formulations detailed above. In Section |4.1[ the 
results are compared against closed-form analytical solutions due to point dipole sources, available in ho¬ 
mogeneous uniaxial media. In Section 4.2 the results are compared against Finite Element Method (FEM) 
results for several selected cases of practical interest in geophysical exploration. These results also examine 
the effect of anisotropy ratios in the surrounding cylindrical layers on the resulting electromagnetic fields. 
Throughout this section, field values are expressed in a phasor form under convention. 
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Figure 5: Relative error distribution with k = 4: (a) Umax = 10^ '^int = 1000, (b) Umax = 20, riint = 1000, 
(c) Umax = 10, flint = 2000, and (d) Umax = 20, flint = 2000. 


4 . 1 . Homogeneous uniaxial media 

Numerical results from the present algorithm are compared to closed-form field expressions due to point 
dipole sources in homogeneous uniaxial media. For a derivation of such analytical solutions, refer to |Ap- 
Ipendix ~K\ The fields are evaluated throughout a square region of 10 cm x 10 cm in the pz-plane. The 


source is a ^-directed Hertzian electric dipole with unit dipole moment and operating frequency of 36 kHz. 
The medium has Cp^h = 16eo [F/m], jUh = 16/io [H/m], ah = 16 [S/m], where eo and /uq denote free-space 
permittivity and permeability values. These horizontal values are fixed whereas different /i^, and ay 
values are considered to yield different anisotropy ratios. 

Fig. - 5d show the relative error between the present algorithm and the closed-form analytical 
solution for different maximum orders Umax employed in the azimuth summation and for various numbers of 
quadrature points Uint employed in the numerical integration. It is assumed that ep^y = eo [F/m], /i^ = /io 


19 



































(a) 




(c) 



Figure 6: Relative error distribution with k = 2: (a) Umax = 10, u-int = 1000, (b) Umax = 20, = 1000, 

(c) Umax = 10, riint = 2000, and (d) rimax = 20, ni„t = 2000. 


[H/m], and cr„ = 1 [S/m], with = k = 4. The relative error is defined as 


relative error^B = 10 log^o 


\Ez,a ~ Ez,n\ 

LsTj 


(98) 


where and indicate analytical and numerical results, respectively. As expected, smaller relative 
errors are obtained for larger number of quadrature points or summation terms. Fig. show the 

relative error distribution for k,^ = k,^ = k, = 2^ under the assumption of Cp^y = 4eo [F/m], fiy = 4/io [H/m], 
and ay = 4 [S/m]. As expected, higher Umax and riint produce smaller relative errors. Fig. 7a -show 


the relative error distribution for = n = v/2, under the assumption of ep^y = 8eo [F/m], fiy = 8/io 

[H/m], and ay = 8 [S/m]. Comparing the cases with k = 4^ k = 2^ k = v/2, we observe that the error 
distribution shows a faster rate of decay along the vertical spatial direction for larger k. Finally, Fig. [8a| 
-[M|show the relative error distribution for k,^ = = k = under the assumption of Cp^y = 16eo [F/m], 

fly = 16/io [H/m], and ay = 16 [S/m], which recovers the isotropic case. In order to further scrutinize the 
effect of Umax and riint ^ the receiver point is next fixed at p — p' = 10 cm, 0 — = 0°, and z — z' = 
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Figure 7: Relative error distribution with k, = \/2: (a) Umax = 10, = 1000, (b) Umax = 20, riint = 1000, 

(c) Umax = 10, flint = 2000, and (d) Umax = 20, flint = 2000. 


cm. Figure - |9d| show the error as fimax and flint vary. When fimax is less than about 15, the relative 
error is not reduced despite the increase in the number of quadrature points. Therefore, fimax should be 
set sufficiently large, otherwise a convergence of results with respect to the number of quadrature points 
would be only relative and the final results would be still inaccurate. For the types of scenarios considered 
here, we have observed that fimax ^ 30 to provide absolute convergence. Furthermore, it is observed that 
convergence is achieved faster for larger anisotropy ratios. This stems from an effective spatial “stretching” 
along the p-direction effected when the anisotropy ratio increases. As can be seen in (15a) and ( |15b| ), 
and k‘^ decrease with an increase in the anisotropy ratio and consequently, higher order modes (with larger 
fi) exhibit a faster decay away at the receiver point. 


4 . 2 . Cylindrically layered seenarios 

In this section, a number of practical cases of interest are considered to illustrate the applicability of the 
algorithm. In all the cases, both the relative permittivity and relative permeability are set to one so 
that = ^p,v = 1 9.nd = 1, whereas the conductivity tensor (and complex permittivity tensor e. 
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Figure 8: Relative error distribution with k = 1: (a) Umax = 10, u-int = 1000, (b) Umax = 20, = 1000, 

(c) Umax = 10, riint = 2000, and (d) rimax = 20, ni„t = 2000. 


see <§) exhibits uniaxial anisotropy where the horizontal resistivity (reciprocal of horizontal conductivity) 
is set to 5 • m and the vertical resistivity is changed, leading to different anisotropy ratios 

Case 1 is depicted in Fig. |lQa| There are three layers, with the first layer representing a metallic mandrel 
with high conductivity, the mid-layer representing a borehole filled with an isotropic fluid, and the outermost 
layer representing the surrounding Earth formation with uniaxial anisotropy. Case 2 is depicted in Fig. [lObi 
where a metallic casing (third layer) is inserted between the borehole and anisotropic formation. Table 
provides the comparison of corresponding results for Case 1 in terms of the square of anisotropy ratios. 
The discrepancy in the magnitude of magnetic fields can be traced to FEM mesh truncation effects: the 
fields obtained by FEM have smaller magnitudes because the Dirichlet boundary condition at the mesh 
boundary moves the ground potential (originally at infinity) closer to the source location. This causes a 
small offset in the results. This is confirmed by Tablewhich shows the relative difference in the computed 
field magnitudes with excellent agreement. Table |5 and [^provide corresponding results for Case 2. 

Case 3 and 4 are depicted in Figs. |lla| and |llb which are the same as Case 2 except for the operating 
frequencies, which for Case 3 is 1 kHz and for Case 4 is 125 kHz. Tables Izlii and provide the corre- 
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sponding results. 


Case 5 and 6 are depicted in Figs. 12a and 12b For Case 5, the borehole is extended to 16" without 
casing. For Case 6, both the transmitter and receiver are positioned inside the formation, which again 
has uniaxial anisotropy. Tables [D and p!3| provide the comparison of corresponding results for Case 5 and 
Case 6 in terms of the anisotropy ratios squared. Tables and show the relative difference in the field 
magnitude for each case. 


The magnetic field magnitude in the y = 0" plane is shown in Figure for Cases 1 and 2. The field 
is plotted in a decibel scale lOlog^^g 1^1 because of the large magnitude variation. In this scale, the small 
differences in magnitude between = 1 and = 10 observed in Tables 3 and 5 are hardly distinguishable. 
On the other hand, these figures clearly show that Case 1 has less confinement of fields within the source 
layer than Case 2 due to the presence of the metallic casing in the latter case, as depicted in Figure pT} 


5. Conclusion 

We provided a robust algorithm for the stable computation of electromagnetic fields in cylindrically 
stratified media with doubly uniaxial anisotropic layers. Range-conditioned integrands, which were originally 
developed for isotropic media, are extended here for uniaxial media. Associated multiplicative factors used 
for the stabilization are expressed as 2x2 matrices in this case. The results show that the formulation is 
indeed stable and have good error controllability. Illustrative scenarios were included to show applicability of 
the proposed algorithm to geophysical exploration problems involving borehole sensors in Earth formations 
with anisotropic responses. 


Acknowledgement 

We thank Halliburton Energy Services for the permission to publish this work, and Dr. Baris Guner for 
kindly compiling some of the comparison data. 

Appendix A. Analytical solution in homogeneous doubly-uniaxial media 

In this appendix, the closed-form analytical expressions used to obtain for electromagnetic fields in 
homogeneous and doubly-uniaxial media are presented. In such media. Maxwell’s equations with 
time-dependence write as 


V X E(r) = icj7i-H(r), (A.l) 

V X H(r) = -icj? • E(r) + J(r). (A.2) 


Permittivity values are complex-valued so as to include conductivities. Eor simplicity, it is assumed that 
the anisotropy ratio for the permeability tensor coincides with that of the complex permittivity tensor, i.e., 
f<ie = To simplify the derivation, we adopt coordinate stretching techniques. Eor a general exposition 
of the coordinate stretching, refer to [all Eg HD]. To begin with, let us consider modified Maxwell’s curl 
equations with stretched coordinates, i.e.. 


V X E(r) = icc;/iH(r), 

V X H(?) = -icj?E(?) + J(?), 


with a modified nabla operator V defined as 


V = 


^ d ^ d ^ d 


(A.3) 

(A.4) 


(A.5) 
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where x, and z are stretched coordinates defined such that 


u ^ u 


I ^du , 

Jo 


(A.6) 


where Su is the corresponding complex stretching variable, and u stands for x, y, or z. In the above, the fields 
and sources are non-Maxwellian but /i and e are scalars, so the medium is isotropic. Using the technique 
in 


(A.3) and (A.4) are rewritten as 

V X 


iuil(detS^ A-H(?), 


V X ( 5 ^ ■ H(r) ) = -iwe 


^detS'^ 


S ■ E(r) + ^detS'^ 


S'- J(r). 


where a dyadic S is defined as 


1 


S = xx[ — ]+yy[ — ]+zz[ — 


1 


1 


Using the relations between the stretched fields to unstretched (Maxwellian) fields. 


E(r) = 5 .E(?), 

H(r)=F'.H(?), 

J(r) = ^det5'^ S' • J(r), 


(A.7) and (A.8) are rearranged as 

V X E(r) = iuj 

V X H(r) = —iuj 


^(dets) U-S 

^fdetS) 


■H(r), 

•E(r)+J(r). 


(A.7) 

(A.8) 

(A.9) 

(A.lOa) 

(A.lOb) 

(A.lOc) 

(A.ll) 

(A.12) 


These two resulting curl equations can be associated with an effective anisotropic medium and represented 
as 


V X E(r) = iujy • H(r), (A.13) 

V X H(r) = -icj? • E(r) + J(r), (A.14) 


which recover the form of the original curl equations ( |A.l[ ) and ( |A.2[ ). Therefore, electromagnetic fields in 
an homogeneous and uniaxial media with E(r) and H(r) can be easily obtained from E(r) and H(r), which 
are solutions in isotropic media with coordinate-stretching, by the transformations expressed in ( |A.lQa ), 
(A.lOb), and (A.lOc). In order to determine the form of the stretching variables relevant to our problem, 
let us examine the effective anisotropic medium obtained above. The constitutive tensors have the form 


where 


^ = 


^detS'^ 
r ^det5'^ 




S'-S' 




= M, 

= eA, 



Sx‘^ 

0 

1 
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Sx 

0 

1 
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0 

Sy 
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= 
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Sx 
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0 
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0 

0 

^x 

Sz - 


(A.15) 

(A.16) 
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(A.17) 























Using two conditions on the stretching varjnbles for uniaxial anisotropy, and using the wavenumber expres¬ 
sion for the modified Maxwell’s equations k = we can set = Sy = and Sz = Consequently,we 

obtain Ji = and e = ^. Next, let us consider the source transformation (A. 10c). If the source is a point 
Hertzian electric dipole like J(r) = —r'), the coordinate stretching should be carefully treated due to 

the presence of the Dirac delta function. The stretched current density is expressed as J(r) = Ila 6{t — r'). 
From the Dirac delta function properties. 


S{r — r') = 


SxSySz 


-5{r-r'), 


(A.18) 


and from (A. 10c), 


a'5{r — r') = ^detS'^ S ■ a (5(r — r') 


0-1 


0 0 


0 s-1 0 
0 0 s7^ 


(5(r — r'). 


(A.19) 


-/ =-i 


Since Sx = Sy = 1 and Sz = k, we have the source transformation a = S -a', and in homogeneous 
isotropic media, the Cartesian field components due to the Hertzian electric dipole source can be written as 


where 


Ex 

ill 6*'=’’ ^ 
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uje 47rr 
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BXY BXZ 
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AZ 

-AY 

-AZ 

0 

AX 
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0 


A = ifc/r — l/r^, 

B = - Siklr^ + 3/f4, 

X = Sx{x' — x) = x' — X, 

Y = Sy{y' -y)=y' -y, 

Z = Sz{z' -z) = k{z' - z), 


k = uj^y.heh/K., 

r = [{x' - xf + {y' - yf + k^{z' - zf] . 


(A.20a) 


(A.20b) 


(A.21a) 


(A.21b) 


(A.21c) 

(A.21d) 

(A.21e) 

(A.21f) 

(A.21g) 

(A.21h) 

(A.21i) 


25 




















=-l 


Applying field transformations, (A.10a) and (A.10b), and source transformation a = S ■ a', we obtain 
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(A.22a) 


(A.22b) 


Finally, applying the coordinate transformations from Cartesian to cylindrical coordinates, we obtain 


where 
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(A.23a) 


(A.23b) 
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(a) (b) 




(c) (d) 

Figure 9: Relative error distribution in terms of various maximum orders Umax and integration points riint 
with the receiver point at p — p' = 10 cm, 0 — 0' = 0°, and z — z' = 1^ cm: (a) /i: = 4, (b) /i: = 2, (c) /i: = \/2, 
and (d) n = 1. 
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Figure 11: (a) Case 3 in the p 2 :-plane and (b) Case 4 in the p 2 :-plane. 
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Table 3: Comparison of magnetic fields in terms of various anisotropy ratios for Case 1. 


Square of 
anisotropy ratio 

Magnetic field [A/m] 
(FEM) 

Magnetic field [A/m] 
(Present algorithm) 

Computing time 
(Present algorithm) 

1 

10.3116 Z98.1899° 

10.5475 Z98.1390° 

10 sec. 

2 

10.2365 Z97.5910° 

10.4723 Z97.5486° 

32 sec. 

5 

10.1565 Z96.9883° 

10.3924 Z96.9612° 

32 sec. 

10 

10.1070 Z96.6327° 

10.3428 Z96.6152° 

31 sec. 


Table 4: Comparison of magnitude difference in magnetic fields for Case 1. 





FEM 

Present algorithm 

between 

= 1 and 

= 2 

0.0751 

0.0752 

between 

= 2 and 

= 5 

0.0800 

0.0799 

between 

= 5 and 

= 10 

0.0495 

0.0496 


Table 5: Comparison of magnetic fields in terms of various anisotropy ratios for Case 2. 


Square of 
anisotropy ratio hz^ 

Magnetic field [A/m] 
(FEM) 

Magnetic field [A/m] 
(Present algorithm) 

Computing time 
(Present algorithm) 

1 

46.6091 Z118.4181° 

46.6303 Z118.4324° 

15 sec. 

2 

46.6099 Z118.4234° 

46.6311 Z118.4381° 

44 sec. 

5 

46.6110 Z118.4283° 

46.6321 Z118.4432° 

44 sec. 

10 

46.6118 Z118.4310° 

46.6329 Z118.4459° 

44 sec. 

Table 6: Comparison of magnitude difference in magnetic fields for Case 2. 




FEM 

Present algorithm 

between hz^ 

= 1 and tzl = 2 

-0.0008 

-0.0008 

between hz^ 

= 2 and = 5 

-0.0011 

-0.0010 

between 

= 5 and = 10 

-0.0008 

-0.0008 
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Table 7: Comparison of magnetic fields in terms of various anisotropy ratios for Case 3. 


Square of 
anisotropy ratio 

Magnetic field [A/m] 
(FEM) 

Magnetic field [A/m] 
(Present algorithm) 

Computing time 
(Present algorithm) 

1 

290.2144 Z127.4332° 

290.2711 Z127.4185° 

7 sec. 

2 

289.9780 Z127.4252° 

290.0564 Z127.4170° 

22 sec. 

5 

289.7678 Z127.4089° 

289.8157 Z127.4175° 

22 sec. 

10 

289.6754 Z127.3988° 

289.6589 Z127.4189° 

22 sec. 

Table 8: Comparison of magnitude difference in magnetic fields for Case 3. 





FEM 

Present algorithm 

between 

= 1 and 

= 2 

0.2364 

0.2147 

between 

= 2 and 

= 5 

0.2102 

0.2407 

between 

= 5 and 

= 10 

0.0924 

0.1568 


Table 9: Comparison of magnetic fields in terms of various anisotropy ratios for Case 4. 


Square of 
anisotropy ratio hz^ 

Magnetic field [A/m] 
(FEM) 

Magnetic field [A/m] 
(Present algorithm) 

Computing time 
(Present algorithm) 

1 

18.7957 Z110.9753° 

18.8074 Z110.9191° 

7 sec. 

2 

18.7959 Z110.9762° 

18.8076 Z110.9200° 

19 sec. 

5 

18.7962 Z110.9770° 

18.8079 Z110.9209° 

19 sec. 

10 

18.7963 Z110.9774° 

18.8080 Z110.9213° 

19 sec. 


Table 10: Comparison of magnitude difference in magnetic fields for Case 4. 




FEM 

Present algorithm 

between hz^ 

= 1 and tzl = 2 

-0.0002 

-0.0002 

between hz^ 

= 2 and = 5 

-0.0003 

-0.0003 

between 

= 5 and = 10 

-0.0001 

-0.0001 
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Figure 12: (a) Case 5 in the p 2 :-plane and (b) Case 6 in the p 2 :-plane. 
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Table 11: Comparison of magnetic fields in terms of various anisotropy ratios for Case 5. 


Square of 
anisotropy ratio 

Magnetic field [A/m] 
(FEM) 

Magnetic field [A/m] 
(Present algorithm) 

Computing time 
(Present algorithm) 

1 

10.7855 Z100.0572° 

10.7857 Z100.0586° 

10 sec. 

2 

10.6723 Z99.4009° 

10.6721 Z99.4024° 

29 sec. 

5 

10.5553 Z98.7180° 

10.5546 Z98.7190° 

29 sec. 

10 

10.4847 Z98.3108° 

10.4839 Z98.3113° 

29 sec. 


Table 12: Comparison of magnitude difference in magnetic fields for Case 5. 





FEM 

Present algorithm 

between 

= 1 and 

= 2 

0.1132 

0.1136 

between 

= 2 and 

= 5 

0.1170 

0.1175 

between 

= 5 and 

= 10 

0.0706 

0.0707 


Table 13: Comparison of magnetic fields in terms of various anisotropy ratios for Case 6. 


Square of 
anisotropy ratio hz^ 

Magnetic field [A/m] 
(FEM) 

Magnetic field [A/m] 
(Present algorithm) 

Computing time 
(Present algorithm) 

1 

8.1259 Z97.0379° 

8.1326 Z97.0341° 

11 sec. 

2 

8.0817 Z96.5262° 

8.0814 Z96.4841° 

32 sec. 

5 

8.0276 Z95.9964° 

8.0271 Z95.9416° 

32 sec. 

10 

7.9939 Z95.6786° 

7.9933 Z95.6240° 

32 sec. 


Table 14: Comparison of magnitude difference in magnetic fields for Case 6. 





FEM 

Present algorithm 

between 

= 1 and hz^ 

= 2 

0.0442 

0.0512 

between tz^ 

= 2 and hz^ 

= 5 

0.0541 

0.0543 

between tz^ 

= 5 and 

= 10 

0.0337 

0.0338 
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Figure 13: Spatial distribution of the magnetic field magnitude on the y = 0" plane at 36 kHz: (a) Case 1 
with = 1, (b) Case 2 with = 1, (c) Case 1 with = 10, and (d) Case 2 with = 10. 


34 























































